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ABSTRACT 

We apply a new, second-order Godunov code, Athena, to studies of the mag- 
netorotational instability (MRI) using unstratified shearing box simulations with 
a uniform net vertical field and a sinusoidally varying zero net vertical field. 
The Athena results agree well with similar studies that used different numerical 
algorithms, including the observation that the turbulent energy decreases with 
increasing resolution in the zero net field model. We conduct analyses to study 
the flow of energy from differential rotation to turbulent fluctuations to thermal- 
ization. A study of the time-correlation between the rates of change of different 
volume-averaged energy components shows that energy injected into turbulent 
fluctuations dissipates on a timescale of Q^ 1 , where Q is the orbital frequency 
of the local domain. Magnetic dissipation dominates over kinetic dissipation, 
although not by as great a factor as the ratio of magnetic to kinetic energy. We 
Fourier-transform the magnetic and kinetic energy evolution equations and, us- 
ing the assumption that the time- averaged energies are constant, determine the 
level of numerical dissipation as a function of length scale and resolution. By 
modeling numerical dissipation as if it were physical in origin, we characterize 
numerical resistivity and viscosity in terms of effective Reynolds and Prandtl 
numbers. The resulting effective magnetic Prandtl number is ~ 2, independent 
of resolution or initial field geometry. MRI simulations with effective Reynolds 
and Prandtl numbers determined by numerical dissipation are not equivalent to 
those where these numbers are set by physical resistivity and viscosity. These 
results serve, then, as a baseline for future shearing box studies where dissipation 
is controlled by the inclusion of explicit viscosity and resistivity. 



Subject headings: Black holes - magnetohydrodynamics - stars: accretion 
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Introduction 



The process of accretion powers a wide range of astrophysical systems, from protostars 
to quasars. In accretion disks, gravitational energy is converted into other forms including 
bulk outflows, heat, a nd ra diation. In the traditional time-stationary thin disk model of 
Shakura fc Syunyaevl (119731 ). the r, component of the stress, r r< p, is proportional to the 
local pressure, r r ^ = aP. The a model assumes that the accretion energy is deposited as 
heat locally and radiated rapidly, providing a relation between disk emissivity and accretion 
rate. While the a model has proven valuable in interpreting many aspects of accretion 
systems, advancing beyond it will require a more detailed understanding of the stress that 
produces angular momentum transport as well as the physical processes involved in the 
subsequent thermalization and radiation of the orbital energy released by those stresses. 

It is now understood that magn etohydrodynamic (MHD) tur bulence generated by the 
magnetorotational instability (MRI) (IBalbus fc Hawleylll99ll . ll998l ) produces significant Maxwell 
stresses, —B r B^j Ani , and Reynolds stresses, pbv r bv^, that account for transport within ac- 
cretion disks. The absence of an analytic theory for MHD turbulence, however, means that 
direct numerical simulations play an essential role in investigating accretion physics. In this 
regard, local simulations, which reduce the problem to the simplest form that can sustain 
MRI-driven turbulence, have proven very useful. The "shearing box" model is a represen- 
tation of a small patch of the disk constructed by boosting to a local co-rotating Cartesian 
frame that ignores geometric cu rvature but re t ains a ll rotational forces. MRI shearing box 
simulations were introd uced by lHawlev etajj (119951) and have bee n extensively used since 



then both wit hout (e.g.. lHawlev et al 



ification (e.g.. iBrandenburg et al 



1996 



19951 : 



Stone et al. 



1996 



Balbus &: Hawlevlll998l) and wit h vertical strat- 



Hirose et all 120061 ) 



Shearing box simulations can investigate several key questions including the functional 
dependence of the stress on disk properties and the turbulent energy flow that leads to 
dissipation as heat. These simulations have made it increasingly clear, for example, that the 
basic a stress parameterization is not only too simplistic, it is actually misleading. Shearing 
boxes have provided ample evi dence that stress is not determined b y pressure, at least in 
the usual manner of the a disk (lHawlev et al.lll995l : ISano et al.l 120041 ) . Early studies showed 
instead that stress is (in some cases) proportional to the magnetic pressure, but the magnetic 



energy is not itself directly determined by the gas and radiation pressure. iBlackman et al. 



( 120081 ) recently reviewed a large number of shearing box results and found that this result 
holds across the full ensemble of simulations with only small differences in the constant of 
proportionality from one run to another. The implications of these results are significant. 
For example, recent local simulations using stratified shearing boxes and radiation transport 
( iBlaes et al.ll2007l ; iKrolik et al.l 120071 ) have found no evidence of the thermal instability long 
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believed to be present in radiation-pressure supported a disks. 



If stress is proportional to magnetic rather than total pressure, what determines the 
magnetic pressure in a disk? Apart from the expectation that the field will remain sub- 
thermal, this remains uncertain. The simplest shearing box simulations using ideal MHD 
have a limited range of significant parameters; this is both a strength and a weakness of 
that model. The magnetic energy in the saturated state could depend upon such factors as 
box size, the amplitude and geometry of the imposed initial ma gnetic field, and the ratio 
of the gas pr e ssure to magnetic pressure (the plasma (3 value). lHawley et al.l (119951 ) and 
Hawley et al.l (119961 ) studied the effect of initial magnetic field topology on the resulting 
stress and found that although the MRI leads to turbulence regardless of the initial field, 
simulations that had an imposed net vertical field produce higher turbulence levels than an 
imposed toroidal field or a s imulation that began with zero net magnetic flux within the 
domain. lHawley et al.l (119951 ) found that the total magnetic energy and the resulting stress 
in the saturated turbulent state was a function of the initial plasma (5 with a uniform vertical 
field, namely that larger j3 (i.e., weaker fields) leads to smaller saturation levels. Other initial 
field configurations do not yield so direct a correlation between background field strength 
and saturation. Many simulations have failed to find any noticeable correlation between 
me an turbulent mag netic energy and the gas pressure. A comprehensive parameter study 
by lSano et al.l (120041 ) observed at best only a very weak gas pressure dependence. 



Since the mean magnetic energy at saturation is presumably a balance between contin- 
ued driving by the MRI and loss due to magnetic dissipation and reconnection, there has 
been interest in going beyond ideal MHD to include explicit physical dissipation in the form 
of kinematic viscosity, u, and Ohmic resistivity, 77. Both of these properties have been shown 
to be imp ortant in dete r minin g; the mean energie s and stresses in MRI t urbulence. Simu- 
l ation s bvlHawlev et al.l (119961). Sano et a l. ( 199 80. IFleming et al.l (120001 ) . ISano fc Inutsuka 



(boOlh . Izieder fc Riidigerl (boOlh . and 



Sano fc Stond (120021 ) have investigated the impact 



of a nonzero rj. The main result of these studies is that increasing the resistivity leads to 
a decrease in turbulence, independent of the initial field configuration. In zero net field 
models, the effect of resistivity on the tu rbulence is larger than one might expect fr om the 
linear MRI relation ([Fleming et al.l 120001 ) . On the other hand, lHawley et al.l (119961 ) found 
that increasing the viscosity increased the magnetic energy in the saturated state. Recent 
work has clarified the situation by demonstrating a dependence of the saturation level on 
both rj and v in terms of the magnetic Prandtl number, P m = u/r]. In particular, the level 
of angular momentum transport increases with increasing P m for simulation s initiated with 



a uniform as well as vanishing mean magnetic field in the vertical direction (IFromang et al. 
20071 iLesur fc Longarettill2007h . 
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Determining the stress levels in MRI turbulence is only one aspect of the problem; 
another is exploring how that turbulence is dissipated into heat. This question has direct 
relevance to phenomenological disk models as well as observations. The a model assumes 
that t he accretion energy is deposited as heat locally and rapidly, and iBalbus fc Papaloizou 
( 119991 ) showed that this property should hold for the energetics of MHD turbulence as well. 
In the simulations, we can determine the rate at which turbulent energy is thermalized and 
the path that energy takes as it moves from the free energy of the shear flow to turbu lence 
and then to heat. Such issues were briefly touched on by iBrandenburg et al.l (119951 ) who 
found that the turbulent magnetic energy was ~ 6 times greater than the perturbed kinetic 
energy, but dissipational heating resulted from roughly equal contributions of magnetic and 
kinetic energy dissipation. This result led t hem to suggest th a t ther e was a net transfer of 
magnetic energy to turbulent kinetic energy. ISano fc Inutsukal (120011 ) studied energy flow in 
the context of MRI channel modes, which are st rong radial streaming motions that result 
from the linear growth of the vertical field MRI (IHawley fc Balbusl Il992l ; IBalbus fc Hawley 
19981 ). Their work included Ohmic resistivity (but not viscosity) and showed that resistive 
heating dominated the thermalization of energy stored in these channel modes. Dissipational 
heating also plays an important role in radiative effec ts and determining di sk structure, both 



of which may be observable properties of disks (e.g., iBeckwith et al.ll2008l ) 



In any study that depends on simulations, there remain factors which cannot be over- 
looked: the effects due to numerics and finite resolution. The majority of the results to- 
date were obtained with numerical codes based on the finite-difference ZEUS algorithm 
dstone fc Normanl Il992al fbh. carried out at relatively low resolution. ZEUS is effectively 
first-order in asymptotic convergence, and in its most widely used form, evolves the internal 
rather than the total energy equation. There have been improvements in both the available 
computational power, which makes higher resolutions and longer evolution times possible, 
and in the algorithms for compressible MHD. In this work, we will reexamine the properties 
of MHD turbulence in the shearing box using a higher-order, Godunov scheme. 



The new code, Athena, (see IStone et al.l 120081 ) represents an improvement over ZEUS 
in several ways includin g true second-order convergence, increased effective resolution (see 
Stone fc Gardiner! 120051 ) . accurate shock capturing, and conservation of total energy. The 
energy-conserving properties of Athena allow us to study energy flow and dissipation within 
the shearing box in greater detail than allowed for by the ZEUS algorithm. The version 
of Athena we use in this paper does not include explicit resistivity or viscosity and instead 
relies on numerical dissipation to thermalize the turbulent energy. Nevertheless, this work 
will serve as a starting point for planned studies of nonideal effects, in cluding the influence of 
P m on the turbulence (IFromang et al.ll2007l ; lLesur fc Longarettill2007l ). As an important part 
of establishing a baseline of simulations, we will characterize the numerical resistivity and 
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viscosity of Athena for the she aring box problem. To do so, we will follow the recent work of 
Fromang Sz Papaloizou (120071 ) who studied the numerical effects of ZEUS on the saturated 



state of MRI shearing box simulations that begin with zero net field. They found that 
the amplitude of the turbulence decreases with increasing resolution and developed several 
useful diagnostics with which to quantify the effective numerical resistivity and viscosity in 
the problem. 

The structure of the paper is as follows. In § [2J we describe the algorithm employed and 
our simulations. In § [31 we reexamine some of the results from previous MRI studies and 
provide a comparison with these studies. In § [U we present the first of two diagnostics used 
to study turbulent energy flow and dissipation. The second of these diagnostics is applied 
in § Finally, we discuss our results and summarize our conclusions in § 



2. Numerical Simulations 

The code used for all of our simulations is Athena, a second-order accurate Godunov 
scheme for solving the equations of ideal MHD in conservative form. The equa tions are solved 
using the dimensionally unsplit corner transport upwind (CTU) metho d of IColellal (119 90') 
couple d with the third-order in space piecewise parabolic method (PPM) of lColella &: Woodward 
( 119841 ) and a constrained transport (C T) algorithm for pr e serving the V • B = co nstraint. 



Det ails of the alg o rithm are described in lGardiner fc Stond (j2005al ). lGardiner fc Stond (120081 ) 
and IStone et al.l (120081 ) . The Athena code has been extensively tested against various hy- 
drodynamic and MHD tests (IStone et al.l 120081 ) . 



We employ the shearing box formalism, in which our computational domain is corotating 
with the fluid flow at some radius in the disk. The domain size is small compared to this 
radius, a llowing us to e xpan d the equations of motion in Cartesian form, as described in 
detail bv lflawlev eTaD Jl995h . 



In the ideal MHD approximation, the evolution of the fluid in the shearing box is 
described by: 



dp 
dt 



dpv 



V ■ (pv) = 0, 
1 



+ V • (pvv - BB) + V(P + -5 2 ) = 2qptt 2 x - 2fl x pv, 



dB 

~dt 



V • (vB - Bv) 



0. 



(1) 
(2) 
(3) 
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— + V • \{E + P + -B 2 )v - B{B ■ v)\ = 2qtfpv ■ x, (4) 

where p is the mass density, pv is the momentum density, B is the magnetic field, P is 
the gas pressure, E is the total energy density, and q is the shear parameter, defined as 
q = — dlnQ / dlnR. Q is the angular velocity of the center of the shearing box. Note that our 
system of units has the magnetic permeability p = 1. We use q = 3/2, appropriate for a 
Keplerian disk. The first source term on the right-hand side of equation (j2J) and the term 
on the right-hand side of equation (@J correspond to tidal forces (gravity and centrifugal) in 
the corotating frame. The second source term in equation (j2J) is the Coriolis force. The total 
energy density is the sum of the thermal, kinetic, and magnetic energy densities 



E 



-pv 2 + -B 2 

r 2 



(5) 



where e is thermal energy density. The equation of state is that of an ideal gas, e = P/(j — l), 
where the adiabatic index is 7 = 5/3 in all simulations. The terms on the right hand sides 
of equations (T5]) and (j3]) are added to the MHD equations in a directionally unsplit manner, 
consistent with the CTU algorithm. Note that we have neglected vertical stratification. 

An important component of shearing box simulations is the sheari ng periodic bo u ndary 
conditions at the x boundaries, which are implem ented as descr i bed in lHawley et al.l (119951 ) 
with a few modifications for Athena. First, as in lHawley et al.l (119951 ). the y momentum is 
adjusted to account for the shear across the x boundaries as fluid moves out one boundary 
and enters at the other. Since Athena evolves the total energy, however, this energy must 
also be adjusted to account for the difference in y momentum across the boundaries. Second, 
following the description in lHawley et al.l (119951 ). quantities are linearly reconstructed in the 
ghost zones from appropriate zones in the physical domain that have been shifted along 
y to account for the shear across the boundary. However, we have found that the precise 
conservation of a quantity depends on how this reconstruction is performed; the fluxes of 
a particular conserved quantity must be reconstructed to conserve the quantity to roundoff 
level. For example, consider the conservation of magnetic flux through the computational 
domain. For the magnetic flux through the box to be conserved to machine precision, 
the line integral of the electromotive force (EMF), £ = —v x B, along the boundaries 
must remain zero. The y and z boundary conditions are periodic, and therefore, the line 
integrated EMFs along these boundaries cancel. This is not the case with the shearing 
periodic boundaries, however. Consider the net B z flux through the grid, which will be 
conserved if £ y = v z B x — v x B z is zero when integrated along both x boundaries. Computing 
the EMF using ghost zone variables v z ,B x ,v x , and B z after reconstruction introduces a 
truncation error, and the B z flux is not conserved. This is avoided if we instead perform 
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the shearing-periodic reconstruction step on S y itself. A similar argument applies to mass 
conservation; one needs to reconstruct the density flux in the shearing boundaries instead of 
the density itselfQ 

In the code, we only perform this EMF/flux reconstruction for S y . We have found that 
conservation of the B z flux is essential, owing to the strong effect on the turbulence due to 
a net vertical field. The perfect conservation of B y is not as important, and as ensuring its 
precise conservation involves a more complex procedure, we allow the B y flux to be conserved 
only to the truncation level. Similarly, the precise conservation of mass has minimal impact 
on the behavior of the turbulence, and we allow the mass to be conserved to truncation 
level. We would like to note, however, that because of this, mass is lost during the MRI 
evolution in our simulations. To quantify the level of mass loss, the total percentage of mass 
lost over 100 orbits of evolution is ~ 2% for our highest resolution simulations (see below 
for a description of our simulations) and ~ 10% for our lowest resolution simulations; we 
observe convergence of mass conservation with resolution. 

Although Athena conserves total energy, the shearing b oundaries do work on the fluid 



and represent a significant energy source. As was shown in lHawley et al.l (119951 ). one can 

12 



integrate the total energy plus gravitational potential energy, E+p<&, where $ = qVL 2 {^—x 2 ) 



over the domain to obtain 

d(E + p$) qVl 



/ (pv x 5v y - B x B y )dydz, (6) 
Jx 



dt L y L z j x 

where L x , L y , and L z are the domain sizes in the x, y, and z directions respectively (see 
below) , 



5v y = Vy + qVtx, 



(7) 



and the integral is calculated over one of the x boundaries. In our simulations, equation 
is satisfied to truncation level with the error coming from the tidal potential source term 
in equation (j4j). It is possible to rewrite this source ter m to guarantee that equation (jSJ) is 
satisfied to roundoff level (see iGardiner fc Stondl2005bl ). but we have found that this makes 
very little difference to how the total energy evolves. 



Gardiner &: Stond (l2005bl ) point out that the source terms in the momentum equa- 



principle, the same argument applies to momentum and energy conservation, but these equations are 
not conserved to machine precision due to the existence of source terms. 
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tion cannot be written in a purely conservative form and that the x and y momenta are 
tightly coupled through these terms. In the hydrodynamic limit the source terms account 
for epicyclic oscillations, and if the epicyclic kinetic energy (see their equation 8) is not con- 
served to machine precision, coupling between long wavelength modes and epicyclic oscilla- 
tion modes can result from truncation error. Over time this coupl ing can artificially i ncrease 



the kinetic energy. To ensure the conservation of epicyclic energy, iGardiner fc Stone! (l2005bl ) 
evolved the angular momentum fluctuations directly rather than the y momentum, casting 
the equations into a form consistent with uniform epicyclic motion. They then employed a 
Crank-Nicholson scheme to evolve the source terms that govern the evolution of the mome- 
tum fluctuations. In MHD, however, oscillatory epicyclic motion is replaced by unstable, 
growing MRI modes. Epicyclic kinetic energy is not conserved and these special techniques 



are not required. Therefore, we use the standard Athena algorithm (e.g.. lStone et al.ll2008l ) 
to evolve the momentum equations. 



As was done in the original shearing box simulations (IHawley et al.lll995l ) our standard 
shearing box has a radial size L x = 1, an azimuthal size L y = 2tt, and a vertical size 
L z = 1. We initialize a velocity flow with v = —qQxy, with q = 3/2, Q = 0.001, and 
—L x /2 < x < L x /2. In an isothermal disk, the sound speed is c s ~ QH where H is 
the scale height. With L z = H, we have c s = L Z Q, and we define the initial pressure as 
P = pVt 2 L 2 z . With p = 1, we have P = 10~ 6 . In this paper, we consider two initial magnetic 
field geometries that are commonly used in shearing box studies. Models labeled NZ (for 
Net Z-field) have an initial uniform vertical magnetic field, B z , and models labeled SZ (for 
Sine Z-field) begin with a sinusoidal distribution of B z and have zero net flux through the 
box. Specifically, we initialize the NZ runs with B = y/2P/Pz, and the SZ runs with 
B = v /2P/(3sm[(2n/L x )x}z. In both cases, we set (3 = 1600. This determines the ratio 
of the vertical box size to the fastest growing linear MRI wavelength as L z /X c ~ 4, where 
A c = 27Ta/16/15|^a|/^, and va is the Alfven speed. To seed the MRI, we introduce random 
adiabatic perturbations to P and p with amplitude 5P/P = 0.01. 

For both of these initial field configurations, we have run a full range of grid resolutions, 
from N x = 16, N y = 32, N z = 16 to the highest resolution used in this study, N x = 128, 
iVy = 256, N z = 128, proceeding by factors of two. All of the simulations were run for a 
total of 100 orbits. 

In addition to the standard shearing box simulations, we have run some additional 
experiments designed to further investigate magnetic and kinetic energy dissipation. First, 
we perform a set of simulations in which we remove the velocity shear and the tidal and 
Coriolis force terms, thus removing the energy source that maintains the turbulence. The 
purpose of these simulations is to investigate energy flow and dissipation in the absence 
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of the shear, which is the driving force for the turbulence. We perform these simulations 
by restarting each of the standard shearing box runs at a time when the shearing periodic 
boundaries are stri ctly periodic. Thes e "periodic points" are given by t n = n L y / qfl L x , with 



n — 0, 1, 2, (see lHawley et al.lll995l ). We choose the restart time to be 40 orbits. We then 



evolve the system to follow the decay of the kinetic and magnetic energies. 

Finally, we run a set of low resolution simulations with varying aspect ratio to examine 
the effect of secondary parasitic modes on the channel solution (see § !3.2p . These simulations 
have the same initial conditions as the net flux simulation with N x = 32, N y = 64, and 
N z = 32 but with varying domain size in the x and y dimensions. The grid cell size (e.g., 
L x /N x in the x direction) in each dimension is kept constant. All simulations are summarized 
in Table [TJ 
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Table 1. MRI Simulations with Athena 



Label 


Initial Field Geometry 


Resolution (N x X N y X N z ) 


Domain (L x x L y x L z ) 


Description 


NZ16 


net flux 


16 x 32 x 16 


1 x 2tt x 1 




NZ32 


net flux 


32 x 64 x 32 


1 x 2tt x 1 




NZ64 


net flux 


64 x 128 x 64 


1 x 2» x 1 




NZ128 


net flux 


128 x 256 x 128 


1 x 2?r x 1 


fiducial run - net flux 


NZD128 


net flux 


128 x 256 x 128 


1 x 2tt x 1 


decaying turbulence 


SZ16 


zero net flux 


16 x 32 x 16 


1 x 2tt x 1 




SZ32 


zero net flux 


32 x 64 x 32 


1 x 2tt x 1 




SZ64 


zero net flux 


64 x 128 x 64 


1 x 2tt x 1 




SZ128 


zero net flux 


128 x 256 x 128 


1 x 2tt x 1 


fiducial run - zero net flux 


SZD128 


zero net flux 


128 x 256 x 128 


1 x 2tt x 1 


decaying turbulence 


NZAR1 


net flux 


16 x 64 x 32 


i X 2tt x 1 


varied aspect ratio 


NZAR2 


net flux 


64 x 64 x 32 


2 X 2tt X 1 


varied aspect ratio 


NZAR3 


net flux 


32 x 32 x 32 


1 X 7T X 1 


varied aspect ratio 


NZAR4 


net flux 


16 x 32 x 32 


ixirxl 


varied aspect ratio 


NZAR5 


net flux 


64 x 32 x 32 


2xixl 


varied aspect ratio 


NZAR6 


net flux 


128 x 32 x 32 


4xixl 


varied aspect ratio 
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3. General Properties of MRI Turbulence 

This work represents the first detailed study of the MRI with Athena, which has an 
algorithm significantly different from that used in ZEUS. To begin, we will reexamine many 
of the shearing box models and the results already documented in the literature. Any 
significant differences between Athena results and those previously published could indicate 
where numerical effects (algorithm, resolution) have an influence. Since Athena is an energy- 
conserving, shock-capturing algorithm it has at least the potential to produce somewhat 
different results. Conversely, agreement between Athena and other codes would support the 
robustness of the shearing box results to date. 

In this section, we describe some of the general properties of MRI turbulence as simu- 
lated with Athena and compare our results with those in the literature. These properties will 
also serve as a starting point for further analysis presented in the following sections. In what 
follows, the highest resolution runs NZ128 and SZ128 will serve as our fiducial simulations 
for each initial field geometry. We study resolution effects for each field geometry using the 
lower resolution simulations. 



3.1. Characteristics of Saturation 

Figures [1] and [2] show the development of the MRI and the subsequent evolution of 
the resulting MHD turbulence for the fiducial NZ128 and SZ128 runs respectively. The 
MRI saturates before orbit 5 and the MHD turbulent state lasts for the remainder of the 
100 orbit simulation. Along with these figures, we list several time- and volume-averaged 
quantities from the fiducial runs in Table [2j The time average is done from orbits 20 to 
100, and the errors are given by one standard deviation over this period. Volume-averaged 
values are indicated by the single-angled bracket notation (e.g., (B 2 )), and time- and volume- 
averaged values are denoted by double-angled brackets (e.g., ((B 2 ))). In both fiducial runs, 
the toroidal field magnetic energy dominates with (Bj/2) > (B 2 £ /2) > (B 2 z /2). Examining 
the components of the kinetic energy and perturbed kinetic energy, which is (p/2)(v 2 + Sv 2 + 
v 2 ) with 5v y given by equation ([7]), we find they are closer to each other in value than are the 
components of the magnetic energy. The relative ordering is similar except that the x kinetic 
energy is larger than the perturbed y kinetic energy, pSv 2 /2, in SZ128. Another feature of 
note is the greater saturation level and fluctuation amplitude of the NZ128 run compared to 
that of SZ128. As in past studies, the Maxwell stress dominates over the Reynolds; the ratio 
of the Maxwell to Reynolds stress oscillates between 1 and 10. Similarly, past studies have 
shown a t ight correlation betwee n Maxwell (and total) stress and the magnetic energy density 
(see, e.g.. lBlackman et al.ll2008l ). Here the ratio of the Maxwell stress to the magnetic energy 
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density is roughly 1/2. These y a lues and the overall observ ations are genera l ly con sistent 
with the results of lHawley et al.l (119951 ). lHawley et al.l (119961 ). and ISano et al.l (120041 ) . 



One major difference from past ZEUS simulations is the evolution of the total (E + p$) 
and thermal (e) energy densities, shown in the lower right plot of Figs. [1] and [2] for the NZ128 
and SZ128 runs respectively. Since we evolve an adiabatic equation of state and there is no 
cooling term in the energy equation, the total energy increases with time at a rate given 
by equation The total energy increases because the free energy of the shearing fluid is 
being thermalized by the turbulence, but the shearing box boundary conditions continuously 
reinforce that shear. The stresses at the radial boundaries therefore constitute a source term. 
Equation ([6]) also explains why the total energy reaches a higher value at the end of the 
simulation in NZ128 compared to SZ128. Since the volume-averaged stress (which is roughly 
equal to the stress at the radial boundaries) is higher in NZ128, the energy injection rate will 
be larger. These plots also show that the thermal energy follows the total energy very closely. 
That i s, the injected energy ends up as thermal energy a short time later (IGardiner fc Stone 
2005bj). We will further study the thermalization of injected energy in §Hand 351 



Does the significant i ncrease in therma l energy affect the turbulence in any way? This 
question was examined by lSano et al.l (120041 ) in an extensive series of simulations. They found 
evidence of a very weak dependence of the time-averaged Maxwell stress on the gas pressure. 
Such an increase is not apparent from a first look at Figs. [TJ and [21 but short timescale 
fluctuations are a dominant feature of these volume-averaged quantities. We examined the 
long term behavior of the Maxwell stress using time-averaging procedures to smooth away 
the fluctuations (which do not appear to change over long timescales). We found marginal 
evidence for a weak dependence of the Maxwell stress on the gas pressure in some, but not 
all, of the data. While it is possible that longer evolution times and a wider exploration of 
parameter space could be useful to address this question further, it is clear the stress has 
barely changed despite an increase in thermal pressure by a factor of order 100 in run NZ128. 
Thus if there is any dependence of the stress on the pressure, it is very weak and does not 
significantly affect the characteristics of local MRI turbulence. 

We study the effect of resolution through a series of lower resolution simulations (see 
Table [TJ. Figure [3] shows the time- and volume-averaged magnetic and perturbed kinetic 
energies as a function of grid resolution for both the net flux and zero net flux initial con- 
ditions. The time average is calculated from orbits 20 to 100; the error bars indicate one 
standard deviation. For the net flux simulati on, there appears to be a slight trend of in- 
creasing energy with resolution, as observed in lHawley et al.l (119951 ). Resolution has a more 
obvious effect on the zero net flux initial condition. The turbulent energies decrease with 
increasing resolution. This resolution effect was previously reported for zero net field initial 
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conditions in other simulations (jFromang &: Papaloizoul 120071 ; iPessah et al.l 120071 ) using dif- 
ferent numerical algorithms. With Athena, the time- and volume- averaged total magnetic 
energy density decreases by roughly a factor of two for each factor of two resolution increase. 
The amplitude of the fluctuations in the total magnetic energy density decreases by roughly 
a factor of two to four for each resolution increment. At all resolutions, the y magnetic 
energy density continues to be the largest, followed by the x energy, and then the z energy. 
As was the case for NZ128, dominates for all net flux simulations, followed by pv^/2, 

and then pv 2 z /2. In the zero net flux simulations, the x kinetic energy density is greater than 
the perturbed y kinetic energy density. These components of the perturbed kinetic energy 
density are close in value, and it is often the case that the x and y components are within one 
standard deviation of each other. The ratio of time- and volume-averaged Maxwell stress to 
time- and volume-averaged magnetic energy density is constant with resolution. The ratio 
of time- and volume-averaged Maxwell stress to time- and volume-averaged Reynolds stress 
has a slight increase with resolution in the net flux simulations and a slight decrease with 
resolution in the zero net flux simulations. However, we point out that the observed trends 
in the ratio of stresses are subject to considerable uncertainty given the large error bars 
calculated for the various quantities. 



3.2. Channel Solution 



One of the interesting aspects of the vertical field MRI in a shearing box is that the 
fastest gro wing mode leads to ax is ymmetric radial stream ing motions, dubbed "channel 
solutions" (IHawley h Baibuslll992l ). iGoodman fc Xul (11994 ) pointed out that for the vertical 
field in an unstratified box, the linear MRI eigenmode is also a nonlinear solution in the 
incompressible limit. They further show that the nonlinear channel solution is itself unstable 
to "parasitic modes." These modes require radial and azimuthal wavelengths larger than 
the vertical wav elength of the channel s olution and will disrupt the channel flow if the box 
is large enough (IBalbus fc Hawleylll998l ). 



In the present simulations, the initial vertical field is sufficiently weak that the fastest 
growing vertical wavelength is less than the radial and azimuthal dimensions of the box, and 
any initial tendency toward the channel solution at the end of the linear growth phase is 
quickly disrupted. However, we find that the large fluctuations in the magnetic energy den- 
sity for NZ128 are a result of recurring channel solutions^ Figure H] shows the azimuthally- 



2 The recurrence of the channel solution presumably results from the fact that the net vertical magnetic 
field can never be destroyed or removed from the domain, given the periodic boundary conditions and the 
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averaged velocities at several times during the amplification and subsequent decay of one 
such fluctuation. The flow organizes itself into a two-channel solution, which becomes more 
well-defined as the magnetic energy incr eases. The channel sol ution is eventually destroyed 
via secondary, parasitic instabilities (see iGoodman &: Xulll994j ). which coincides with a de- 
crease in magnetic energy. The same channel solution appears during other instances of 
large magnetic energy fluctuation in NZ128 and does not appear in SZ128. Furthermore, 
the recurring channel flows appear in the lower resolution net magnetic flux simulations. 
As observed previously, the channel solution an d large magnetic energ y fluctuations are a 
property of simulations with a uniform B z field ( ISano fc Inutsukall200ll ). 



Since the channel solution is subject to parasitic modes that depend on the available 
wavelengths that can fit in the box, we expect that this behavior is influenced by the domain 
aspect ratio employed. To verify this, we have run several low resolution simulations (labelled 
NZAR1 — NZAR6, see § [2]) using different aspect ratios. We found that for large enough 
L x , the int e rmitte nt channel modes no longer occur; this behavior was also observed by 
Bo do et al.l (120081 ). The prominence of intermittent channel flows is a consequence of the 
restrictions introduced by the domain size. However, we use this property in § HI where the 
large fluctuations in turbulent energy created by the channel solutions provide a clear marker 
of energy injection by the boundaries. We can then track the subsequent thermalization of 
that energy. 



3.3. Energy Power Spectra 

The nature of MRI-driven MHD turbulence can be characterized in part by the power 
spectrum of kinetic and magnetic energies. To obtain such power spectra, w e do a full 3D 



Fourie r transform on the simulation data employing the procedures outlined in lHawley et al. 



( 119951 ) to account for the shearing-periodic boundaries. Briefly, the shearing periodic bound- 
ary conditions in the x direction allow the domain to be strictly periodic in the x direction 
only at certain times, called periodic points t n (described in §2]). To perform a standard fast 
Fourier transform (FFT) at some time t that is not equal to t n , we transform the data into a 
frame where the x boundaries are strictly periodic. We then calculate the FFT in this frame 
and remap to the original frame. 

The turbulent magnetic, kinetic, and perturbed kinetic energy densities in Fourier space 
are defined as 



strict conservation of z magnetic flux. 
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\B x (k)\ 2 + \B y (k)\ 2 + \B z (k)\'< 



i|^(fc)| 2 = \ [lv^Ofc)l 2 + IVM,(*0I 2 + IV/^Wi 2 



(9) 
(10) 



where / means the Fourier transform of / defined by 



/(*) 




/(*)e 



ik-X-,3 



d 6 x. 



Note that for the kinetic energies, we include the density along with the velocity when 
calculating the Fourier transform, resulting in the appearance of J~p in the above equations. 
To obtain these quantities as a function of length scale and to improve statistics, we average 
our data over shells of constant k = \k\. For further improvement of statistics, we average 
each of these terms over 161 frames (i.e., from orbit 20 to 100 in increments of 0.5 orbits). 

Figure [5] shows the power spectra of these energy densities for the net flux and zero net 
flux runs. The figure shows resolution effects as different lines in each plot. In all cases, the 
largest scales account for most of the energy. The general shape of the energy powe r spectra 
agrees with previous studies (e.g. jHawley et al.lll995t iFromang fc Papaloizoull2007l ). For the 
net flux simulations, the magnetic energy dominates over the kinetic and perturbed kinetic 
energies at all scales, independent of resolution. As the resolution is increased, the power 
spectra extend to higher k, but the general shape remains constant. At some values for k, the 
uncertainty in energy (not plotted), represented by one temporal standard deviation around 
the mean, is large enough to overlap with other energy components, making it difficult to 
conclusively say which energy dominates at these particular scales. 

We calculated a power law index in Fourier space for each energy density and at each 
resolution. This slope was determined by a linear fit to the energy densities in log space from 
fcL/(27r) = 1 to the maximum scale for the given resolution. There is some uncertainty in 
this measurement because the power spectra are not strictly linear in log space (see Fig. [5]). 
In NZ128, the energy density is proportional to [kL/(2ir)} n with n m —4 for every energy 
density. This index is approximately constant with resolution, but there is evidence that n 
becomes more negative at higher resolutions. In determining an error in the value of n, we 
found that this error is often dominant. Thus, such a resolution dependence is somewhat 
tentative. 
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There is a noticeable resolution dependence in the zero net flux simulations. First, as 
resolution is increased, the magnetic energy density decreases at all scales. This effect was 
discussed in $ 13.11 the power spectra are consistent with the power spectrum analysis of 



Fromang fc Papaloizoul (120071 ). The same resolution dependence is observed for the per- 
turbed kinetic energy density. The magnetic energy density at small k decreases faster with 
resolution than does the perturbed kinetic energy density. The total kinetic energy density 
(i.e., including shear) remains constant with resolution, which simply results from the fact 
that the shear velocity, which dominates the kinetic energy, is constant with resolution. The 
uncertainty in each energy component appears to be smaller than in the net flux simulations. 
However, there are still some values of k at which the calculated errors overlap. 

We calculated a power law index in Fourier space for each energy density and resolution 
for the zero net flux simulations. The procedure we used was the same as for the net flux 
simulations. For the kinetic and perturbed kinetic energy densities, we found that n lies 
between -3.5 and -4, whereas for the magnetic energy density, n lies between -3 and -3.5. 
There does not appear to be any resolution dependence in n for the magnetic energy density, 
but there is a tentative decrease in n (similar to the net flux case) with increasing resolution 
for the kinetic and perturbed kinetic energy densities. 
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Table 2. Saturation Characteristics 





Quantity 


NZ128 


SZ128 






((-B X By))/P 


0.216 ± 0.116 


(6.55 ± 1.15) x 10" 


-3 




{(pV X 5Vy))/P 


0.028 ± 0.019 


(1.91 ± 0.76) x 10" 


-3 




«B 2 /2»/Po 
{(Bl/2))/P 


0.488 ± 0.262 


0.014 ± 0.003 






0.071 ± 0.027 


(2.01 ± 0.38) x 10" 


-3 




({Bl/2))/P 


0.388 ± 0.231 


0.011 ± 0.002 






((Bl/2))/P 


0.029 ± 0.011 


(7.98 ± 1.57) x 10" 


-4 




((p5vi/2))/P 


0.145 ± 0.060 


(7.69 ± 1.81) x 10" 


-3 




{(pvl/2))/P 


0.046 ± 0.024 


(3.73 ± 1.27) x 10" 


-3 




{(p5v y /2))/P 


0.078 ± 0.035 


(2.68 ± 0.60) x 10" 


-3 




((pvl/2))/P 


0.021 ± 0.011 


(1.28 ± 0.21) x 10" 


-3 


«- 


B X By))/((pV X 5Vy)) 


7.60 ± 6.47 


3.43 ± 1.49 




« 


-B x B y ))/{(B*/2)) 


0.443 ± 0.336 


0.462 ± 0.116 
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4. Energy Fluctuations 

Athena evolves the equation for total energy, the volume-average of which will change 
only due to the Maxwell and Reynold stresses at the radial boundaries (equation (jSJ)). As 
was discussed in §|3l the individual volume-averaged magnetic and kinetic energies are highly 
variable throughout the evolution as energy is continuously transferred between magnetic, 
kinetic and thermal components. We can study these energy flow processes by tracking 
the energy injected at the boundaries as it is subsequently thermalized in the turbulence. 
For this purpose, the existence of the recurring channel solution in the net magnetic field 
simulation is very useful; the sudden increase in stress provides a clear injection of energy 
that can be traced using several diagnostics. Having developed these diagnostics we can then 
apply them to the zero net magnetic flux simulations. Finally, to gain additional insight into 
dissipation in the turbulence, we conduct an experiment in which the shear flow and gravity 
terms have been removed. 



4.1. Sustained Turbulence 

The total energy density, including the gravitational potential energy density, is defined 

as 

E tot = E + p$ = e + l -pv 2 + l -B 2 + p$ (12) 

where $ is given in £j2j Averaging equation (jT2l) over the entire domain, taking the time 
derivative, and rearranging the terms, we obtain, 

f — E- m — K — M — G. (13) 

where E in = d(E tot )/dt is the energy injection rate due to stress at the boundaries (see 
equation T = d(e)/dt is the rate of change of thermal energy density, K = d(^pv 2 ) /dt 
is the rate of change of the kinetic energy density, M = d(^B 2 ) /dt is the rate of change of 
the magnetic energy density, and G = d(p&)/dt is the time derivative of the tidal potential 
energy density. The brackets indicate a volume-average over the simulation domain. G is the 
change in a fluid element's gravitational energy as it moves within the domain. We expect 
the contribution of the tidal potential term to be insignificant, an expectation borne out by 
direct computation. We will ignore this term in most of the subsequent discussion. The 
stress terms at the radial boundaries are generally positive, which means energy is being 
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injected into the box via the work done by this stress (E- m > 0). The remaining terms in 
equation (fT3l) can be either positive or negative. 

The lower right plot in Figure [T] shows that the thermal energy density closely follows 
the total energy density, but with a short time delay. This can be better seen in Figure El 
which shows the individual terms from equation (|T3|) for a 20 orbit period in the NZ128 
simulation. There is a clear time delay of less than one orbit between significant changes in 
the energy injection rate and the thermal energy derivative, suggesting a comparab le delay 
befor e the injected energy is ther malized, a property noted in lSano fc Inutsukal (120011 ) as well 



as m iGardiner & Stone! (l2005bl ). These features in the energy derivatives result from the 
creation and destruction of channel flows. During this time interval, the magnetic and kinetic 
energies are also changing. By examining the maxima in the thermal energy derivative and 
the corresponding features in the kinetic and magnetic energy derivatives, it appears that 
the magnetic energy dissipation dominates the thermalization process. 

It is useful to define a temporal correlation function for the various energy components 
by writing 



< N-\L\-1 



c 



N - L 



/] A i+ \ L \Bi 



i=0 



i N - 1 



if L < 



AB 



N 

i=0 
N-L-l 



(14) 



N- L 



A ^ 



+L 



i=0 



JV-1 



if L > 



i=0 



where A and B are two time-series datasets N elements in length. The quantity L is the 
number of elements over which to shift A and B with respect to each other to calculate the 
correlation coefficient. We apply equation f[T4"|) to the energy rates by setting A = T, and 
B = K, M, or E in . This allows us to correlate the energy injection rate and the change in 
kinetic and magnetic energies against the change in thermal energy over certain timescales. 
Since T > 0, if the correlation between T and K (or M) is negative, then kinetic energy (or 
magnetic energy) must be decreasing, and a strong negative correlation would suggest that 
kinetic energy (or magnetic energy) is being thermalized. 

Figure [7] is the correlation function for B = E- m , K, and M calculated over orbits 20 to 
100. The x-axis is the correlation timescale in units of orbits. We only look at correlation 
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times of < 1 orbit as the degree to which the thermal energy evolution follows that of the 
total energy (see Figures [1] and [2]) indicate that thermalization happens over that timescale. 
To examine the correlation function on longer timescales would be misleading since peaks 
in the function would suggest a correlation between two events that are not causally related 
(e.g., the injection of energy for one channel event being correlated with the thermalization 
of energy for another channel event). The left plot of the figure shows that E- m is strongly 
correlated with T on a timescale of At ~ -0.2 orbits. This correlation is exactly what we 
observed in Fig. O The energy injected by the stress at the boundaries ends up as heat 
less than one orbit later. The negative sign on this value of At simply means that the 
injection happens before the thermalization. In the right plot, both K and M are negatively 
correlated with T suggesting that magnetic and kinetic energy are being thermalized. The 
stronger magnetic correlation further suggests that magnetic dissipation contributes more to 
thermalization than kinetic dissipation. The positive correlation between K and M against 
T at negative At values is a result of the magnetic and kinetic energies increasing along 
with the energy injection into the box. That is, the stress at the boundaries increases the 
magnetic and kinetic energies which are dissipated a short time later. 

An interesting feature is evident in Fig. [7J the negative peak in the magnetic and kinetic 
correlation functions occur for At slightly greater than zero. Similarly, in Fig. [6] one can see 
that peaks in the magnetic and kinetic energy derivatives are offset with respect to the energy 
injection and thermalization peaks. For example, the maximum rate for magnetic energy 
loss occurs after the maximum rate for thermal energy gain. Of course, these are plots of 
the time derivative of the energy, so a peak simply indicates where the second derivative is 
zero. The magnetic energy is both losing energy to dissipation while gaining energy from the 
shear at the boundaries. When the energy injection rate peaks decline, the thermalization 
rate is still growing and the magnetic energy rate also peaks and begins to decline. Similarly, 
the slope of the magnetic energy loss rate will change sign after the thermalization rate has 
peaked and when the energy injection rate is no longer itself in decline. 

As a test, we performed this correlation analysis on the lower resolution net-flux simu- 
lations and find that energy injection precedes thermalization by ~ 0.2 orbits, independent 
of resolution. Furthermore, magnetic dissipation dominates over kinetic dissipation for all 
net flux simulations. 

The analysis so far has only examined the rate of change in the energy terms, not specif- 
ically how they change. For example, does a "dip" in M correspond to direct thermalization 
of magnetic energy, or is there a transfer of energy from magnetic to kinetic? To examine the 
energy flow in more detail, we focus on orbits 50 to 52 in NZ128, for which we ran the NZ128 
simulation at high temporal resolution. This high time resolution allows us to resolve short 
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timescale features, but also generates many large data files. Therefore, we restrict this part 
of the analysis to the two orbit period mentioned above. Consider the evolution equation 
for the volume-averaged kinetic energy given by 



K = — < V 



v Q/w 2 + ^B 2 + P + p$J - B(v ■ B) 
+ ( ( P + ^B 2 j V-vJ-(B-(B- V«)> - G - Q k , (15) 

where Qk is the volume-averaged (numerical) kinetic energy dissipation rate. The evolution 
equation for the volume-averaged magnetic energy is given by 

^ = ~ ( v ■ (l B2v ) ) ~ {l B2v ■ v ) + (B ■ (B ■ Vv)) " ^ m (16) 

where Q m is the volume-averaged (numerical) magnetic energy dissipation rate. We have 
calculated each term in these equations over the two orbit period and find that the dominant 
terms are —(V • (\pv 2 )v), (V • [B(v ■ B)]), (B ■ (B ■ Vv)), Qk, and Q m . Q k and Q m are what 
remain after calculating all other terms in the energy equations at a particular instant in time. 
Calculating the volume-averages of th e first two terms yie lds the radial boundary Reynolds 



and Maxwell stresses in Equation (jBJ) ( IHawley et al.lll995l ). namely the energy injection rate 



by the shearing periodic boundaries. The third of the dominant terms is the transfer rate 
of kinetic to magnetic energy via field line stretching. Figure [8] plots the time-history of 
this term (pink line) along with T (black line), the energy injection rate E- in (blue line), and 
— Qk and — Q m (green and red lines, respectively). As energy is injected into the grid, a 
significant fraction of this energy is transferred to the magnetic field via field line stretching, 
presumably through the shear flow. Thermalization follows 0.2 orbits later and is marked 
by increases in the absolute value of Qk and Q m , with \Q m \ > \Qk\- The ratio of kinetic to 
magnetic dissipation is approximately constant in time over this period, with Qk/Q m ~ 0.6. 
This suggests that the details of the thermalization do not vary with intermittent increases 
in E- m that occur when the fluid experiences a channel flow. 

As discussed, the recurring channel modes in the net flux simulations create distin- 
guishable points of energy injection that make it straightforward to follow the subsequent 
thermalization. Such modes do not exist in the zero net flux simulations, which makes the 
identification of specific correlations slightly more difficult. The situation is further com- 
plicated by the overall reduced levels of the turbulence which causes the time derivative of 
the thermal energy to be dominated by very high frequency oscillations due to propagating 
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spiral density waves ([Gardiner Stone! l2005bl ). We have determined that these waves are 



created by compressibility and have very little effect on the dissipational heating within the 
box. To remove their dominance in the energy derivatives, we rebin the time data using a 
"neighborhood" averaging procedure in which the rebinned data points are calculated from 
averages of a specified number of original data points. We then apply equation (fl~4"j) between 
E- m and T; the result is shown in Fig. [9j The correlation curve has several narrow peaks, 
which result from residual effects of the rebinning process. The curve has a broader peak 
near At ~ -0.2 orbits, which agrees with the same curves for NZ128 (Fig. [7]). The correlation 
function for SZ128 is not as sharply peaked as that for NZ128, which is most likely a result 
of the lower amplitude variability in the rebinned SZ128 data. Applying this analysis to the 
lower resolution zero net flux simulations, we find that the correlation function always has a 
broad peak at At ~ -0.2 orbits. Thus, as was the case in the net flux simulations, the energy 
injection/thermalization timescale is independent of resolution. 

Finally, we note that the saturated state of SZ128 is too complex to obtain correlations 
between M, K, and T, such as was done for NZ128. In the net flux simulations, the 
recurring channel modes lead to the build up and thermalization of magnetic energy. The 
creation and thermalization of magnetic energy are events that are well-separated in time, 
making it easy to study the flow of energy between various components. In the zero net 
flux simulations, however, the average properties of the turbulence remain more constant in 
time. We will further investigate the dissipation of magnetic and kinetic energy for the zero 
net flux geometry in §4.21 and §5.11 



4.2. Decaying Turbulence 



As noted by lHawley et al.l (119951 ). the MHD turbulence decays without differential ro- 
tation to sustain the MRI. We make use of this to observe how rapidly thermalization occurs 
when there is no further input of energy. This analysis should provide some additional in- 
sight into the thermalization process for each field geometry. We remove the net shear flow 
and the Coriolis and tidal forces from a state taken from the sustained MRI turbulence in 
the fiducial models. These runs are labeled "NZD" and "SZD" in Tableland are described 
in more detail in £j2j Figure [10] shows the subsequent magnetic and kinetic energy decay for 
both runs. In the figure, the kinetic and magnetic energies have been normalized to their 
values at the starting time of t = 40 orbits. 

In NZD128, the ratio of total magnetic to kinetic energy at t = 40 orbits is 3.4. The 
figure shows that the magnetic energy decays more rapidly than the kinetic energy at early 
times, losing almost half its initial value within 0.2 orbits. In SZD 128, the ratio of total 
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magnetic to kinetic energy at t = 40 orbits is 1.4. The kinetic energy shows high frequency 
oscillations about an average value that decays in time. These oscillations are due to the same 
compressive, spiral waves that exist in the sustained turbulence simulations. The magnetic 
energy is unaffected by these waves. The average decay of kinetic energy, calculated from 
smoothing away the oscillations, is also shown in the figure. Both the kinetic and magnetic 
energies decay quickly over time. Again, almost half the magnetic energy is lost within 0.2 
orbits. The high frequency oscillations also decay in amplitude over time. As was the case 
in NZD128, the magnetic dissipation rate is initially faster than that for the kinetic energy. 
After about one orbit, the decay rates become comparable. 

Finally, we checked the contributions from the various terms in equations f|T5|) and (fT6j) . 
In both NZD128 and SZD128, there is some transfer from magnetic to kinetic energy during 
the decay. However, the transfer rate is small compared to the decay rate of the magnetic 
energy and is such that the numerical dissipation of magnetic energy dominates over that of 
kinetic energy. 



Transfer Functions 



In their investigation of convergence of zero net flux shearing box simulations, iFromang fc Papaloizou 



(120071 ) carried out an analysis based on the evolution of magnetic energy in Fourier space. 
This analysis shows how magnetic energy is created, transferred from one scale to another, 
and finally lost due to numerical dissipation. Their study used the ZEUS code and as- 
sumed an isothermal equation of state. Here we repeat and expand upon their analysis to 
understand dissipation as a function of length scale in Athena. 



We note several differences between our work and that of IFromang &: Papaloizou! (120071 ) . 
First, they focus on magnetic energy evolution and did not provide a comparable calculation 
for the kinetic energy. Second, recognizing that the y direction is dominated by the largest 
scales, they restricted their analysis to axisymmetric modes, namely k y = 0. Finally, as they 
were primarily interested in how poloidal field could be regenerated as part of a dynamo 
process, a portion of their analysis concentrated on t he poloidal components rather than the 
full magnetic field. We have chosen to extend the IFromang fc Papaloizou! (120071 ) analysis 
more generally to include a kinetic energy density evolution, nonaxisymmetric effects, and 
the effects of a nonzero toroidal field. 



Following IFromang fc Papaloizou! (120071 ). we decompose the velocity field of the flow 
into the mean flow, V s h, and the turbulent velocity, v t , via 
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v = V sh + v 



(17) 



The mean flow is denned as 



V sh = V sh y 



y 



LyL Z 




v y (x,y,z)dydz. 



Turning next to the induction equation, we substitute equation ffTTj) for the velocity, 
take the Fourier transform, and dot the result with the complex conjugate of B(k), which 
is defined by equation ( TTTT) with f = B. All Fourier transforms are done via equation ( TTTi) 
using a standard FFT and replacing / with the appropriate quantity. The data is mapped 
into a frame in which the x boundaries are periodic and then remapped into the original 
frame after performing the FFT. 



The result of this calculation is an equation describing the magnetic energy density 
evolution in Fourier space, 



id\B(k)\ 2 

2 dt 



A + S + T bb + T divv + T bv + D 



(19) 



where 



A = -Re 



B*(k) 




dy 



(20) 



S = +Re 



Bl{k) 




Bx ^ e - ikx d 3 x 
ox 



(21) 



T, 



hi, 



-Re 



B*{k) 




v t ■ V)Be- lk - x d 3 x 



(22) 



- divt; 



-Re 



B*(k) 




(V • v t )Be~ ik - x d 3 x 



(23) 



T bv = +Re 



B*(k) 




B ■ V)v t e~ ik - X d z x 



(24) 
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The -D mag term has no analytic expr ession; it is simply what is le ft over and accounts for 



numerical losses of magnetic energy (IFromang fc Papaloizoul 120071 ). In the present simula- 



tions, there is no physical resistivity The other terms have the following meanings: A is the 
transfer of magnetic energy between scales by the shear flow, S is the creation of magnetic 
energy from this shear flow, Tbb is the advection of magnetic energy between scales by the 
turbulent velocity field, T$ nv results from the turbulent compressibility, and Tf, v describes the 
creation of magnetic field by the turbulent velocity fluctuations. In each case, Re signifies 
the real part of the transform. 

One can follow a similar procedure using the momentum equation to determine the 
evolution of the kinetic energy density in Fourier space. As described previously, we include 
the density in our Fourier transforms. Consider the time derivative of y/pv given by 

Note that here, for simplicity, we do not decompose the velocity into mean and turbulent 
components. Using a combination of the continuity and momentum equations, this equation 
can be written as 



dyTpv 
dt 



-v -Vv- -V(P + -B 2 ) + -(B • VB) - 2fl x v + 2qVt 2 xx 
p 2 p 



v 



+ — [-p(V ■v)-vVp] 



(26) 



If we take the Fourier transform of this equation and dot the result with the complex con- 
jugate of 



y/pv{k) 




y/p(x)v(x)e- lkx d 3 x, 



(27) 



we arrive at 



l d\y/pv(k) 

2 dt 



T vv ~\~ T com p -\- T v b -\- Tp ress -\- T cov -\- -\- D 



kin; 



(28) 



where 



l -\^pv{k)\ 2 = \ v 



\yfpv x {k)\ 2 + \^fpv y {k)\ 2 + \^/pv z (k) 



(29) 
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T = -Be 




(30) 



T = —Be 

± comp 1 LO 




^(V -v)e-* k x d 3 



x 



(31) 



T vb = +Be 




— {B ■ V)Be~ ik ' x d 3 x 



(32) 



T = —Be 

- 1 press 1 Lc 




-L V ( P+ l B 2 )e- lkx d 3 x 
VP * 



(33) 



cor 1 ^ 




(2fl x y/pv)e- ik - x d 3 x 



(34) 



7^ = +i?e 




2 ifc-a^g. 



(35) 



and -Didn accounts for the dissipation of kinetic energy. Again, this dissipation is numerical 
as we have not included an explicit viscosity term in our equations. Equation ( 1281) describes 
the evolution of the kinetic energy density in Fourier space. T vv is a term that describes the 
transfer of kinetic energy between scales by the velocity field (both the mean and turbu- 
lent velocity), T comp results from turbulent compressibility, T vb describes how kinetic energy 
changes from magnetic tension, T press represents the effect of both gas and magnetic pressure 
on the kinetic energy, and is the effect of the tidal potential on the kinetic energy. Note 
that T cor is analytically equal to zero, and it is not included in any of the following analysis 
or discussion. 

In the saturated state of the MRI, the magnetic and kinetic energy densities should be 
in a steady state on average (although they do show strong fluctuations over short periods 
of time). If we consider the time-averages of equations (fT9l and ff28l) . then we can set the 
left hand sides to zero. We then rewrite these equations as 



T vv ~\~ Tcomp T v f) -\- -7~p ress -|- Tfj) -\- Dki n 0, 



(36) 



A + S + Tbi, + Tdi vv + Tb v + -D mag — 0, 



(37) 
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where each of these terms is now a time-average. Here we average over 161 snapshots from 
orbit 20 to 100 in increments of 0.5 orbits. Each of these terms is a function of k x , k y , and 
k?., and in what follows we av erage the terms on shells of constant k — \k\ as was done in 



Fromang &: Papaloizoul (120071 ) j 3 l Note that unlike the averaging described in that paper, we 



include k y in the calculation of k. 



5.1. Zero Net Magnetic Flux 

5.1.1. Fiducial Run 

In this section, we focus on the Fourier transfer functions for the fiducial zero net 
magnetic flux simulation. Figure [TT] plots the magnetic transfer functions defined in equa- 
tions (|20p - (|24p as a function of length scale for SZ128, and Fig. [15] plots the kinetic transfer 
functions defined in equations (!30l) - (l35l) . The dashed lines correspond to plus or minus one 
standard deviation around the mean value of the time average. Most of the transfer func- 
tions show large variation at small k values which may be due to poor statistics at small k 
and relatively large time variability. Because the transfer functions approach zero rapidly, 
we plot the ranges 1 < kL/(2ir) < 20 and 20 < kL/(27r) < 64 in the same figure, but with 
different y scalings. 



The shear term S is positive at all scales, as observed in iFromang &: Papaloizoul ( 120071 ) 
meaning that B y is creat ed by the shear flow at al l scale s. A is small at all scales, supporting 



the assumption made in IFromang &: Papaloizoul (120071 ) that A w 0. T bv is primarily neg- 



ative at the largest scales, although there are large fluctuations, and becomes positive for 
kL/(2n) > 35. The turbulent velocity fluctuations seem to be creating magnetic energy at 
the smallest scales, but at larger scales, the magnetic field appears to lose energy via this 
interaction with the turbulence. T bb is negative for k smaller than kL/(2ir) ~ 20, meaning 
that the turbulence is transferring magnetic energy away from these scales. Although this 
analysis doesn't determine the direction of this cascade, at the largest scale (i.e., the box 
size) the energy can only cascade to smaller scales. In terms of absolute value, S and T bb are 
dominant on the largest scales, while on small scales, T bb > T bv > S > 0. 

It is difficult to say anything conclusive about the kinetic transfer functions on the 
largest scales as they are subject to considerable uncertainty, although T vb < appears 
reasonably well constrained at these scales. At smaller scales, the two dominant terms are 
T vv and T vb , with T vb > T vv > 0; kinetic energy is being transferred to these scales by the 



3 In our analysis, the average over shells of constant k was done before the temporal average. 
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turbulence, and being created by magnetic field. 



Equations (I3"6"j) and (13"T|) have been set to zero from the assumption that the magnetic 
and kinetic energies are in a time-averaged steady state. The dissipation terms D mag and 
-Dkm are simply what is left over after the other transfer functions have been computed. The 
top plots in Fig. [13] are the kinetic and magnetic dissipation and the ratio D^/ D mag as a 
function of k for 20 < fcL/(27r) < 64; the scatter at small k is large and there is considerable 
uncertainty in the dissipation values. At small scales, magnetic dissipation dominates kinetic 
dissipation by a factor of roughly three. The kinetic and magnetic dissipation rate increases 
in magnitude towards larger scales. 



Following iFromang &: Papaloizoul (120071 ). we can determine an effective resistivity and 
viscosity as a function of length scale by assuming that the numerical effects behave as if they 
were physical resistivity and viscosity. For example, with a constant Ohmic resistivity, the 
induction equation would have an additional term proportional to V 2 B, with the constant 
of proportionality being the resistivity. If we take the Fourier transform of this term and dot 
it with the complex conjugate of B(k), the real part is 



T„ = +Re 



B*(fc). 




d d x 



-k 2 \B(k)[ 



(38) 



We can then define an effective resistivity as a function of k by 



Ves(k) 



-^mag \k) 



(39) 



Similarly, a constant kinematic shear viscosity would add a term proportional to 
v /p[V 2 f+|V(V-f)] to equation (126|) . with the constant of proportionality being the viscosity. 
Note that we only consider shear viscosity here for simplicity. We take the Fourier transform 
of the viscous term, dot it with the complex conjugate of equation ( 12~T1) . and take the real 
part. The result is 



T v = +Re 




^Tp[V 2 v + ^V(V • v)] e - tkx d 3 x 



(40) 



This equation can be made simpler by realizing that the second term of the integrand, 
related to the divergence of v, is negligible. We can also assume that the density is relatively 
constant, and arrive at 
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T u = -k 2 \^pSv(k)\ 2 . (41) 

We have substituted the perturbed velocity here because it is the only velocity that can lead 
to numerical dissipation of kinetic energy. That is, a pure shear flow will not encounter any 
numerical viscosity, and we can subtract off this flow. We define an effective viscosity by 

-<*> - w- (42) 

We can also characterize the effective resistivity and viscosity in terms of a Reynolds 
number, 

Re cS {k) = **L, (43) 

and magnetic Reynolds number, 

Rm cS (k) = -^L (44) 

where we have used the initial isothermal sound speed, c Q = 0.001, as a characteristic velocity, 
and H = L z is a characteristic length. These numbers quantify the numerical dissipation 
coefficients in a dimensionless manner. 

Finally, we define an effective Prandtl number by 

p m Mk) = ^§ ^ 

The effective viscosity and resistivity as well as the effective Prandtl number are shown 
in the bottom plots of Fig. [131 The viscosity and resistivity are fairly constant at large k. The 
effective Reynolds numbers are on the order of Re e s ~ 12000, and Rm e ff ~ 20000 at large 
k. The Pra ndtl number is also relatively flat at these scales, and P m ,es ~ 1-6. This result 



agrees with iFromang fc Papaloizoul (120071 ). where P m:C s > 1 for ZEUS. While the numerical 
dissipation of Athena is not physical, the "flatness" of u e g and rj e g suggests a resemblance to 
physical dissipation at small scales. 

Finally, note that although the Prandtl number is greater than unity, the magnetic 
dissipation dominates over kinetic dissipation. Evidently, T v is larger than T v because there 
is more magnetic energy than kinetic energy at a given scale. In particular, 
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£ = jfo 1 ' . (46) 

T ^ (^^(fc)! 2 

Since there is more magnetic energy than perturbed kinetic energy at a given scale, magnetic 
dissipation dominates. 



5.1.2. Resolution Effects 

To gauge the effect of resolution on these various quantities, we perform the same 
analysis on the lower resolution runs, SZ16, SZ32, and SZ64. We focus, in particular, on 
the small scales (i.e., large k) where our quantities are statistically more well-determined. 
Figure [H] shows u e s, T] e s, P mte s, and the ratio of -Dkin to -D mag as a function of x resolution, 
N x . The data points are calculated by averaging the quantity of interest over k in the regions 
of £;-space where the error on the quantity is less than its mean valueQ The displayed error 
bars are the propagation of the errors from the temporal statistics. At these large values 
of k, v e ff, Ves, P m ,es, and D kin / D mag are relatively flat, varying by a factor of at most 2. 
Consequently, these averages should be representative at small scales. 

The numerical viscosity and resistivity decrease as a function of resolution. The dashed 
lines in the two upper panels of the figure show the line z/ eff ,?7 eff oc N~ 2 . The viscosity and 
resistivity decrease slower than this with increasing N x ; we measured v e ff,Vef{ oc N~ 1-6 . The 
figure also shows that both the effective Prandtl number and the ratio of kinetic to magnetic 
dissipation are constant with resolution to within the error bars. 



5.1.3. Comparison with Previous Results 



Fromang fc Papaloizoul (120071 ) were interested in the transfer function for the poloidal 



field, as the regeneration of this field is key to a self-sustaining dynamo. They found that 
the magnetic dissipation of ZEUS for the poloidal magnetic field departs from the physical 
dissipation model at small k and could even be a nonphysical "positive" dissipation. We 
repeat the same analysis as performed in that paper, but with SZ128, for comparison. First, 
we examine the magnetic dissipation for the full 3D Fourier analysis described above. Second, 
we do the same procedure but setting B y = to focus on the effect of only including 



4 There are some quantities for which the error is never less than the mean. In these cases, we average 
over regions where the mean is greater than 80% of the error. 
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poloidal field. Finally, we perform the procedure with B y = and in the plane k y = 
i. e., axisymmetry) . Thes e simp lifications allow us to reproduce the poloidal field analysis 



of IFromang fc Papaloizoul (120071 ). 



The results are shown in Fig. [151 The left two plots correspond to the Fourier analysis in 
which only B y = is assumed. The right plots assume B y = and k y = 0. The black lines in 
the bottom two plots correspond to the magnetic dissipation for the full 3D Fourier analysis 
with B y ^ and k y ^ 0. It is apparent that when B y = is assumed in the calculations, 
the magnetic dissipation becomes positive at large scales. However, when B y is included, 
the magnetic dissipation remains negative. Whether or not k y = is assumed seems to 
make very little difference, supporting the notion that small k y dominates. Since Athena 
and ZEUS both find positive -D mag at small k, it is unlikely that this effect can be attributed 
to algorithmic limitations specific to ZEUS. Since -D mag is not a derived quantity but simply 
what remains after all the transfer functions are calculated, it seems likely that the positive 
-Dmag values for the poloidal field analysis are due to incomplete statistics at large scales, or 
other inadequacies of the analysis when applied solely to the poloidal field. At small k, the 
standard deviations of the quantities (dashed lines) are considerable. The standard deviation 
on D mag when B y ^ is significantly larger than when one sets B y = 0. This reflects the 
large variability of (By/2) compared to the other components of magnetic energy (see e.g., 
Fig. [2]). At any given time, -D mag can be positive; the assumption of time-stationarity does 
not hold at any point in time. But when the data are time- averaged, -D mag < 0. 

Finally, we compare the numerical magnetic Reyolds number calculated with equa- 
tion (144)) but with the B y = and k v = assumptions. Fo r SZ12 8, we find that Rm c s ~ 



11000, and for SZ64, Rm e $ ~ 3500. IFromang fc Papaloizoul ((20071) find Rm efi ~ 30000 for 



their N x = 128 run, and Rm e R ~ 10000 for their N x = 64 run; both of their calculated 
effective Reynolds numbers are larger than those calculated for Athena. This result seems 
to suggest that ZEUS is actually less dissipative than Athena. However, there are several 
points to consider. First, numerical dissipation is a nonlinear function of resolution, sharply 
increasing as the number of zones per wavelength decreases (high wavenumbers) . The effec- 
tive Reynolds n umber is obtained by measuring d issipation at the high k end of the spectrum. 



As reported by IShen. Stone. &: Gardinerl (120061 ) Athena appears to have higher dissipation 



than ZEUS for poorly resolved waves, as evidenced by the ability of Athena to avoid the 
aliasing errors seen with ZEUS for hydrodynamic shearing box waves. They further point 
out that for wavelengths larger than 16 grid points Athena is less dissipative. Further, 2D 
simulations of decaying turbulence have demonstrated that when saturation amplitude is 
reached, the decay time is longer in Athena than in ZEUS, consistent with Athena having a 



higher effective resolution (IStone &: Gardinerl 120051 ). In the present context, we find that the 



time- and volume-averaged total stresses in our simulations are larger than those calculated 
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in the simulations of iFromang fc Papaloizoul (120071 ). Stronger turbulence leads to larger ki- 



netic and magnetic turbulent fluctuations, which in turn enhances dissipation via grid-scale 
effects. Finally, we reemphasize that assuming B y = may have a significant impact on the 
measurement of effective magnetic dissipation via this analysis. 



5.2. Net Magnetic Flux 

5.2.1. Fiducial Run 



We perform the same transfer function analysis on the fiducial net magnetic flux run, 
NZ128. The various transfer function terms as a function of k are shown in Figs. [ToT - TiTl As 
was the case in the zero net flux simulation, S is positive at all scales and dominates at small 
k; A is relatively small throughout. T bv and T bb are negative at large scales and positive at 
small scales, with T bb > for kL/(2ii) > 5, and T bv > for kL/(27i) > 20. At small scales, 
> Tb v > S > 0. Of the kinetic terms, T vv and T v b dominate with T v b > T vv > 0. These 
results are in general agreement with SZf 28, except that the magnitude of the various terms 
is larger for NZ128 than for SZf28, and and Tb v become positive at smaller k values 
compared to SZf 28. 

As before, we calculate the kinetic and magnetic dissipation as well as effective values 
for the viscosity and resistivity. Figure [TBI shows these quantities for NZf28 at the smallest 
scales. As was the case for SZf28, the mean magnetic dissipation dominates over kinetic 
dissipation by a factor of roughly three at these scales. Note, however, the large error bars 
associated with these plots, which encompass values of -Dkin/Anag > f- Again, the error 
bars are the temporal standard deviation of the transfer functions. Since NZf 28 has a larger 
temporal variability, larger error bars are expected. The mean value for D^ in / D mag is on the 
order of 0.6-0.7, which is consistent with the analysis in § 14. ll in which we found Q^/Q m ~ 0.6. 

The effective viscosity and resistivity show the same basic result as in the SZf28 case. 
v e E, f?efr, and P m , e s change by a factor of order unity at large k. The effective Reynolds 
numbers are on the order of Re c g ~ 4000, and Rm e s ~ 8000 at large k. P mjC fr has a mean 
value of ~ f.9. Again, there is considerable uncertainty in these values due to the large 
amplitude fluctuations in the turbulence. The error bars encompass values of P m ,eS less than 
unity. As a result, it is more difficult to conclusively say that the dissipation behaves the 
same way in NZf28 as in SZf28. However, in an average sense, the two simulations agree 
well qualitatively. 
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5.2.2. Resolution Effects 

We can again look at the effect of resolution on these various dissipation quantities. 
Figure [H shows this effect for the net flux simulations (NZ16, NZ32, NZ64, and NZ128). 
The procedure by which to average over k is the same as described in § 15.1.21 The displayed 
error bars are the propagation of the errors from the temporal statistics. At these large 
values of k, u e s, rjeg, P m , e s, and D^ in / D mag are relatively flat, varying by a factor of at most 
2. 

The numerical viscosity and resistivity decrease as a function of resolution. The dashed 
lines in the two upper panels of the figure show the line u e s,VeS oc Nf 2 . The viscosity and 
resistivity decrease slower than this with increasing N x ; we measured v e ft,Veft oc Nf 1 ' 3 . The 
figure shows that the effective Prandtl number is constant with resolution to within the error 
bars. There appears to be a slight increase in D^/ D mag with resolution, but this trend is 
not definitive given the large uncertainties on the data. 

One might expect i^s and i] e ff to decrease with increasing resolution since these terms 
arise from truncation error. Linear wave advection test problems with Ath ena have shown 



that the truncation error converges at second order (e.g.. lStone et al.ll2008l ). On this basis, 
one would expect u eS) r] e{i oc N~ 2 . We find a shallower decrease with N x , but MRI turbulence 
is a fully nonlinear system, and one should not necessarily expect the same convergence 
behavior as in a linear system. 



6. Discussion and Conclusions 

We have carried out a series of local, unstratified shearing box simulations with the re- 
cently developed Athena code to study the characteristics of MRI driven turbulence. Athena 
uses a second-order, conservative, compressive MHD algorithm, which is significantly differ- 
ent from the algorithms employed in many of the previous MRI studies. In our work, we 
have run several standard models for comparison with previous work, and characterized the 
numerical dissipation of the Athena code for the shearing box problem. Furthermore, we 
have exploited the energy conservation property of Athena to carry out a study of energy 
flow within MRI-driven turbulence. 

To compare with previous numerical results, we have investigated the effects of different 
initial field geometries (uniform or sinusoidal B z ), varying domain aspect ratio, and numerical 
resolution. In all of our simulations, the MRI is initiated and sustained over many orbits. 
The time- and volume-averaged properties of the resulting turbulent flow, such as stress 
levels and magnetic and kinetic energies, are consistent with previous results. As in previous 
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work, we find that boxes containing net vertical field saturate at higher amplitudes compared 
to those without net fields. The total stress is proportional to the magnetic pressure with 
a constant of proportionality ~ 0.5, but is independent of the gas pressure. In the net 
field simulation, the gas pressure increases by a factor of 100, due to thermalization of the 
turbulence, without affecting the stress. The consistency of these results with past work 
indicate that these properties do not result from details of the employed algorithm. 

Fourier analysis of the turbulence shows that the largest scales in the box dominate the 
energetics. In the presence of a net field, the amplitude of the spatial power spectra is largely 
independent of resolution on the largest scales. This is not true for the zero net flux simula- 
tions however. For those simulations, the amplitude decreases as resolution increases, which 
is consistent with the overall resolution behavior. For net field simulations, the averaged tur- 
bulent magnetic and kinetic energies increase slightly with resolution, whereas for the zero 
net field simulations, the energies decrease with increasing resolution roughly in proportion 
to the grid zone size. This apparent lack of convergence for the zero net fi eld shearing box 



simulations was previously demonstrated by iFromang fc Papaloizoul (120071 ) using the ZEUS 
code. 

The net field simulation shows intermittent channel flows which cause temporary in- 
creases in stress through amplific ation of large-scale MRI modes. The parasitic modes de- 



scribed by I Goodman fc Xul (119941 ) destroy the channel flow within about one orbit of time, 
but the rapid increase in stress produces a subsequent increase in thermal energy. The pres- 
ence of these discrete channel flow events is a consequence of the box size — larger boxes 
do not experience them — but we use their presence to study the subsequent energy flow 
following a rapid increase in stress. 

Because Athena evolves the total energy equation, magnetic and kinetic energy losses 
due to numerical grid-scale effects are added to the internal energy. This makes Athena well 
suited to examining the turbulent energy flow and subsequent dissipation. The recurring 
channel flows in the net flux model provide a sudden injection of energy into the box by 
increasing the stress operating on the shearing boundaries of the box. The injected energy 
appears as heat after ~ 0.2 orbits. This corresponds to a timescale Q~ l , which equals L z /c s 
where c s is the initial soundspeed. This timescale determines the amplitude of the Alfven 
speed, va, and its fundamental MRI wavelength, Amri; L z /c s ~ Amri/^a- The timescale is 
thus on the order of the eddy turnover time, indicating that dissipational heating is a local 
process and that energy is not carried over large distances before it is thermalized. 

In the fiducial zero net magnetic flux simulation, SZ128, there are no recurring channel 
modes, making it more difficult to trace the flow of injected energy. The analysis is further 
complicated by the presence of compressive waves that dominate the time derivative of the 
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thermal energy, T. These waves are also present in the net field simulations, but their 
amplitude is smaller relative to the larger turbulent kinetic energy found with a net field. 
A detailed examination of the components of the internal energy equation indicate that 
the compressive waves do not appear to contribute significantly to irreversible heating. By 
averaging T for the zero net flux simulation, we find a correlation of T with E^ n on the same 
timescale of ~ 0.2 orbits. 

In the net field simulation, the dissipation of magnetic energy is larger than that for 
the kinetic energy, not unexpected as the ratio of the average magnetic to perturbed kinetic 
energy is ~ 3.4. But the ratio of the magnetic to kinetic dissipation rate is roughly constant 
at ~ 1.7. The fact that the ratio of dissipation rates does not equal the ratio of energies 
may result from a couple of possibilities. First, there could be a ne t trans fer of magnetic 
to perturbed kinetic energy as was suggested in iBrandenburg et al. Second, the 



difference in the ratios could arise from the effective Prandtl number being larger than one. 
In particular, if oc i / e ffSv 2 /2 and Q m oc r] c sB 2 /2, then (5 2 /5u 2 )(Qk/Qm) ~ P m ,es- With 
the above values for the energy and dissipation ratios, we find (B 2 / ' 8v 2 ){Q^/Q m ) ~ 2, which 
is consistent with the determination of P m , e s from the Fourier analysis (see discussion below). 
The agreement between the two separate calculations of P m ,eS may be coincidental, but it is 
suggestive of Qk oc v eS 5v 2 /2 and Q m oc r] e sB 2 /2. 

The turbulence is sustained by the continued action of the MRI in extracting energy 
from the differential rotation. This can be removed from the simulations allowing us to study 
the decay of the turbulence in detail (simulations NZD128 and SZD128). Figure [TD] shows 
that magnetic losses dominate over kinetic losses during this decay. In both simulations 
nearly 50% of the magnetic energy and 20% of the kinetic energy has been dissipated after 
0.2 orbits. By one orbit into the decay, most of the magnetic and kinetic energy has been 
lost. Although these decay timescales arise in a turbulent flow that lacks power input from 
the MRI, the results are consistent with the conclusion that turbulent energy dissipation 
occurs on a rapid timescale of order Q^ 1 . 



Fromang fe Papaloizoul (120071 ) used a detailed Fourier analysis (§ EJ) to study magnetic 
energy flow and thermalization as a function of length scale in the shearing box. In this 
analysis, the individual terms in the evolution equation for the magnetic energy are exam- 
ined in Fourier space. Averaging over time and assuming that the magnetic energy is in 
a statistical steady state, one sets the sum of these terms equal to a remainder, which is 



5 § 14. II shows that there is in fact a net transfer of kinetic to magnetic energy. However, this kinetic energy 
includes the shear flow, and thus, this result tells us nothing of the energy transfer between magnetic and 
perturbed kinetic energy. 
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credited to numerical effects. These numerical losses can then be modeled as an effective 
resistivity (and viscosity for the kinetic energy), allowing one to characterize the numerical 
dissipation in the simulation. 

We repeated their analysis with Athena and extended it to the kinetic energy. The 
dominant effect at large scales is the generation of magnetic field by the background shear. 
This energy is transferred to other scales by the turbulence. Net positive field creation by 
the turbulent flow and energy gains by the transfer between scales only happens at small 
wavelengths. This point of transition from loss to gain happens at smaller scales for the 
zero net field simulation compared to the net field model. Magnetic dissipation dominates 
over kinetic dissipation at small scales (i.e., kL/{2it) > 20). Modeling these as an effective 
resistivity rj and viscosity v shows that rj and v drop with increasing resolution with a power 
that lies between first- and second-order in grid resolution. The effective Prandtl number, 
on the other hand, is nearly constant as a function of resolution with a value between ~ 1.5 
and 2. 



Fromang fc Papaloizoul (120071 ) observed what they described as "negative" resistivity 
in an analysis restricted to the poloidal field alone. In repeating their exact analysis with 
Athena, we also observed such an "anti-dissipation" at large scales. This indicates that this 
effect is not associated with a numerical algorithm limitation associated with ZEUS. More 
likely, it arises from the statistical uncertainty at large scales and from the failure of the 
assumptions that go into the definition of the dissipation term. We note that the inclusion 
of the toroidal field B y in the analysis shows net dissipation at all scales, although again the 
statistical variation is large at large scales. 

In conclusion, w hat do these results imply for shearing box simulations and the MRI? 



s imp 

First, as observed by lFromang fc Papaloizoul (j2007l ). the scales over which turbulent energy 



generation occurs are not well-separated from those where there is significant dissipation; 
the MRI operates over a wide range of scales. The MRI grows at a rate ~ kv a for all 
k less than Q/va- At large scales, a weak field will grow more slowly than the timescale 
over which energy is transferred between scales, between magnetic and kinetic forms, and 
ultimately thermalized. If a field is chopped up by reconnection, it may be reduced to 
small scales where the MRI no longer operates. In the presence of a net field, there will 
always be a significant driving term at the scales set by that imposed field. In the absence 
of such a field, however, the outcome will be determined by the complex interplay of loss 
due to dissipation and amplification by the MRI. In the numerical simulations with zero 
net field, increasing the reso lution causes an overall decrease in the saturation energies. 



Fromang &: Papaloizoul (120071 ) attribute this to higher resolution enabling the MRI to operate 



at intermediate scales which facilitates the transfer of energy to small scales and promotes 
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reconnection and dissipation. What is perhaps surprising is that resolving the MRI at 
these scales leads to greater field dissipation than would otherwise be accomplished by the 
numerical losses that would occur if those scales were underresolved. Because the same 
effect is observed with both Athena and ZEUS, it seems likely that this ability of the MRI 
to transfer energy away from the largest scales in the shearing box and to increase the total 
dissipation is a physical rather than numerical effect. 



In related work, iFromang et al.l (120071 ) and iLesur fc Longarettil (120071 ) studied the effect 
of varying the physical (not numerical) magnetic Prandtl number, P m , on the turbulence. 



They found that the saturation amplitudes were increased with increased P m . IFromang et al. 
( 120071 ) found evidence that there exists a critical P m > 1 below which zero net field sim- 
ulations would die out rather than achieve a steady turbulent state. Our results in this 
investigation show that this Prandtl number dependence is a distinct effect from the ob- 
served dependence of the turbulence on resolution. We find the numerical P m to be largely 
independent of resolution in Athena. Taken together, however, the dependence on physical 
P m and the dependence on resolution point to the importance of small and intermediate scale 
magnetic dissipation and reconnection to establishing saturation amplitudes in MRI-driven 
turbulence. 



As discussed by IFromang &: Papaloizoul (120071 ). numerical dissipation can deviate sig- 
nificantly from physical dissipation. In § 15.1.11 we showed that r) e s and u e g are relatively flat 
at small scales, suggesting a resemblance to physical dissipation. However, consider the nu- 
merical Reynolds number as calculated from equation (|43p for our zero net flux simulations. 
For N x = 128, we found Re e s ~ 1200 0, and P m , e ff <y 1.6 for all of zero net flux simulations. 



From the parameter space studies of IFromang et al.l (120071 ). these values for the Reynolds 
and Prandtl numbers correspond to marginal MRI turublence; that is, they lie very close to 
the critical line between sustained and decaying turbulence. For N x = 64, Re e s ~ 4100, and 
the Reynolds number is even smaller for the lower resolutions. These values are well within 
the decaying turbulence regime, but we find active MRI turbulence in all of our simulations. 
These results show that the effective Reynolds and Prandtl numbers of Athena as measured 
at large wavenumbers does not apply at smaller k values where there are many grid zones 
per wavelength. Thus, the Reynolds numbers and Prandtl numbers that we calculate should 
be taken as a measure of the effective numerical dissipation of the code and not equated 
to a flow with the same Reynolds and Prandtl number as determined by a simple physical 
resistivity and viscosity. 

This result highlights an uncertainty associated with any MRI simulation that depends 
only on numerical rather than physical dissipation. It is apparent that the numerical Prandtl 
number can play an important role in determining the ratio of magnetic to kinetic dissipation. 
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More speculatively, the Prandtl number may also play a role in the timescale over which 
thermalization occurs. In the present study, we found that both the thermalization timescale 
and the effective numerical Prandtl number were largely independent of resolution. However, 
the turbulent energy thermalization timescales and properties we measure may be subject 
to change when explicit dissipation is included. It will be a very important next step in this 
work to include physical dissipation and verify these results. 

This work is only the first step in applying Athena to the problem of the energetics of 
MRI turbulence. The present study provides a calibration of the numerical dissipation, which 
will be important in future studies that include explicit resistivity and viscosity. Furthermore, 
the unstratified shearing box has the virtue of simplicity and allows a detailed study of MRI 
turbulence without too many confounding factors, but it also may prove too limited for 
predictive application to accretion flows. The inclusion of vertical stratification and radiative 
cooling are both straightforward extensions to the present study. The detailed diagnostics 
developed and applied in this study should prove valuable in this planned work. 

We thank Jim Stone, Steve Balbus, and Sebastien Fromang for useful discussions and 
suggestions regarding this work. We also thank the anonymous referee whose comments and 
suggestions improved this paper. This material is based upon work supported by NASA 
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Fig. 1. — Volume- averaged energy densities and stresses normalized to the initial gas pressure 
versus time for the NZ128 simulation. In the upper two plots, the black line is the total energy 
density, the green line is the component of the energy density in the x direction, the red line 
is the y direction component, and the blue line is the z direction component. The upper 
left plot shows the volume-averaged magnetic energy density, the upper right plot shows the 
perturbed kinetic energy density (i.e., with the shear subtracted off of v y ), and the lower left 
plot is the volume-averaged total stress (black), Maxwell stress (pink), and Reynolds stress 
(blue). The lower right plot is the total energy density, including gravitational energy (solid 
line), and the thermal energy density (dashed line). The y axes have the same range for all 
plots except for the total/thermal energy density plot. 
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Fig. 2. — Volume- averaged energy densities and stresses normalized to the initial gas pressure 
versus time for the SZ128 simulation. In the upper two plots, the black line is the total energy 
density, the green line is the component of the energy density in the x direction, the red line 
is the y direction component, and the blue line is the z direction component. The upper 
left plot shows the volume-averaged magnetic energy density, the upper right plot shows the 
perturbed kinetic energy density (i.e., with the shear subtracted off of v y ), and the lower left 
plot is the volume-averaged total stress (black), Maxwell stress (pink), and Reynolds stress 
(blue). The lower right plot is the total energy density, including gravitational energy (solid 
line), and the thermal energy density (dashed line). The y axes have the same range for all 
plots except for the total/thermal energy density plot. 
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Fig. 3. — Time- and volume-averaged energy densities normalized to the initial pressure for 
various resolutions. The two upper plots correspond to the net flux simulations, and the 
two lower plots correspond to the zero net flux simulations. The left plots are the averaged 
magnetic energy densities, and the right plots are the averaged perturbed kinetic energy 
densities (i.e., with shear subtracted off of v y ). In all plots, the black symbols are the total 
energy density, the green symbols are the x component of the energy density, the red symbols 
are the y component, and the blue symbols are the z component. The time averages are done 
from orbit 20 to 100, and the error bars indicate one standard deviation over this period. 
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Fig. 4. — The development and destruction of a channel flow during the NZ128 simulation. 
The upper left plot shows a fluctuation in the volume-averaged magnetic energy density from 
t = 80 orbits to t = 85 orbits. The remaining plots show the y-averaged perturbed y 
velocity (colors) and v x and v z (vectors). The upper right plot occurs at t = 82.5 orbits, 
the lower left plot occurs at t = 83 orbits, and the lower right plot occurs at t = 84 orbits. 
These times are indicated on the upper left plot by the arrows. At t — 82.5 orbits, one can 
see the development of a two-channel flow, in which one channel has v x < and 5v y < 0, 
and the other channel has v x > and 5v y > 0. At t = 83 orbits, this channel flow is 
even more developed as the perturbations to the y velocity have become even stronger and 
v x dominates over v z everywhere. The development of this channel flow coincides with an 
increase in volume-averaged magnetic energy density. By t = 84 orbits, the channel flow has 
been destroyed, coinciding with the decrease in magnetic energy density. 
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Fig. 5. — Spatial power spectra of various energy densities in the saturated state of the 
standard net flux (left panels) and zero net flux simulations (right panels). The spectra were 
obtained via an average over 161 frames in the saturated state and an average over shells 
of constant modulus |fc|. In each column, the first plot shows magnetic energy density, the 
second shows kinetic energy density, and the third shows perturbed kinetic energy density 
(as defined in the text). The effect of resolution is shown in each individual plot; the dotted 
line corresponds to the resolution with N x = 16, the dot-dashed line corresponds to N x = 32, 
the dashed line corresponds to N x = 64, and the solid line corresponds to N x = 128. All 
energy densities have been normalized to the initial gas pressure and are plotted against a 
dimensionless wave number (L is the length of the smallest dimension of the box). 
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Fig. 6. — Time derivative of various volume-averaged energy densities for a 20 orbit period 
in the NZ128 simulation. The time derivative of the energy densities have been multiplied by 
an orbital time over the initial gas pressure. The dark blue line is the energy injection rate, 
E- m , the black line is the thermal energy density derivative, T, the green line is the kinetic 
energy density derivative, K, and the red line is the magnetic energy density derivative, M. 
The dotted line indicates zero. 
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Fig. 7. — Correlation coefficients calculated over orbits 20 to 100 in the NZ128 simulation. 
The plot on the left was calculated by correlating the energy injection rate, E in , against the 
thermal energy time derivative, T. The x-axis is the correlation length in time, and the 
y-axis is the coefficient multiplied by an orbital period over the initial gas pressure. The plot 
on the right was calculated by correlating the magnetic energy derivative, M, (solid line) 
and kinetic energy derivative, K, (dashed line) against the thermal energy derivative. The 
dotted line indicates Cab = 0. Note that the two plots have different y scales. 
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Fig. 8. — Various terms in the volume-averaged magnetic, kinetic, and thermal energy 
density evolution equations over a two orbit period of NZ128. The energy terms are T 
(black), E in (blue), -Qk (green), -Q m (red), and the volume-averaged transfer rate from 
kinetic to magnetic energy (pink). All of these terms are defined in the text and have been 
multiplied by an orbital period over the initial gas pressure. 
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Fig. 9. — Correlation coefficient calculated over the saturated state of the SZ128 simulation. 
The coefficient was calculated by correlating the energy injection rate against the thermal 
energy density derivative. The x-axis is the correlation length in time, and the |/-axis is the 
coefficient multiplied by an orbital period over the initial gas pressure. The narrow peaks in 
the curve correspond to residual effects from rebinning the energy derivatives (described in 
the text). The broader peak in the correlation function occurs at At ~ -0.2 orbits. 
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Fig. 10. — Volume-averaged magnetic and kinetic energy densities in the first 1.5 orbits 
of NZD128 (left) and SZD128 (right). In both plots, the upper curves correspond to the 
kinetic energy density and the lower curves correspond to the magnetic energy density. In 
SZD128, high frequency oscillations appear in the kinetic energy evolution. To smooth away 
these oscillations, a moving window average was applied to the kinetic energy density. The 
unsmoothed kinetic energy is shown by the dotted line, while the smoothed kinetic energy is 
the solid line. The magnetic energy density in the SZD128 plot has also been smoothed for 
consistency. Both the kinetic and magnetic energy densities have been normalized to their 
respective (unsmoothed) values at t — 40 orbits. 



-51 - 



S ond A (x10" M ) 



1.0 
0.5 
0.0 

-0.5 
-1.0 



S 









0.0030 
0.0015 
0.0000 

-0.0015 
-0.0030 




1.0 
0.5 

0.0 

-0.5 
-1.0 



0.0030 
0.0015 

0.0000 

-0.0015 
-0.0030 




k L/(2ti) 



Fig. 11. — Magnetic Fourier transfer functions versus a dimensionless wave number (L is 
the length of the smallest dimension in the box) for SZ128. Each plot is displayed in two 
components; the left part shows the data for 1 < kL/ (27r) < 20, and the right part shows the 
data for 20 < hL/(2ir) < 64 by changing the x and y axis scaling. In all plots, the solid line is 
the average value for the transfer function. This average was obtained over 161 frames in the 
saturated state and shells of constant |fc|. The upper (lower) dashed line that matches color 
with the solid line correspond to the transfer function plus (minus) one temporal standard 
deviation. From top to bottom, the plots show S (red) and A (black), T bb , T divv , and T bv . 
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Fig. 12. — Kinetic Fourier transfer functions versus a dimensionless wave number (L is 
the length of the smallest dimension in the box) for SZ128. Each plot is displayed in two 
components; the left part shows the data for 1 < kL/ (27r) < 20, and the right part shows the 
data for 20 < kL/ (2tt) < 64 by changing the x and y axis scaling. In all plots, the solid line is 
the average value for the transfer function. This average was obtained over 161 frames in the 
saturated state and shells of constant |fc|. The upper (lower) dashed line that matches color 
with the solid line correspond to the transfer function plus (minus) one temporal standard 



deviation. From top to bottom, the plots show T v 



T vb , T press (red) and T comp (black), and 
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Fig. 13. — Numerical dissipation quantities plotted against a dimensionless wave number (L 
is the length of the smallest dimension in the box). These plots correspond to data from 
SZ128. The upper left plot shows the dissipation rate of kinetic energy (green) and magnetic 
energy (red) in Fourier space. The upper right plot shows the ratio of these two dissipation 
rates. The lower left plot shows the effective numerical viscosity (green) and resistivity (red). 
The lower right plot shows the ratio of the viscosity to resistivity (i.e., the effective Prandtl 
number). In all plots, the solid line is the average value for the quantity of interest. For 
.Dkin and D mag , this average was obtained from averaging over shells of constant |fe| and 
over 161 frames in the saturated state. The averaged viscosity and resistivity values were 
calculated as described in the text. The upper and lower dashed lines correspond to the 
error propagated from one temporal standard deviation. 
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Fig. 14. — Averaged dissipation related quantities as a function of grid resolution. These 
plots correspond to data from the zero net flux simulations, SZ16, SZ32, SZ64, and SZ128. 
The upper left plot shows the effective viscosity versus x resolution. The dashed line shows 
u c g oc N~ 2 . The upper right plot shows the effective resistivity versus x resolution. Again, 
the dashed line shows i] c g oc N~ 2 . The lower left plot shows the ratio of kinetic to magnetic 
dissipation versus x resolution. The lower right plot shows the effective Prandtl number 
versus x resolution. For each resolution, the data point was obtained from averaging the 
quantity as a function of k over values of k where the error in this quantity is not much larger 
than the mean value. The error bars represent the propagated errors from the temporal 
statistics. 
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Fig. 15. — Magnetic dissipation rate from the SZ128 simulation for three versions of the 
transfer function analysis. The upper left plot and the red lines in the lower left plot corre- 
spond to the analysis in which B y = was assumed. The upper right plot and the red lines 
in the lower right plot correspond to the analysis in which both B y = and k y = were 
assumed. The black lines in the lower plots result from relaxing both of these assumptions. 
The solid lines in the upper plots correspond to -D mag whereas the dashed lines correspond 
to rjTjj with 77 = 10~ 7 chosen to provide a reasonable match to D mag at large k. The dashed 
lines in the lower plots correspond to one standard deviation above and below the quantity 
represented by the solid line of the same color. A horizontal line at zero is shown in all plots 
as the blue dotted line. Note the difference in ?/-axis scale between the upper and lower 
plots. 
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Fig. 16. — Magnetic Fourier transfer functions versus a dimensionless wave number (L is 
the length of the smallest dimension in the box) for NZ128. Each plot is displayed in two 
components; the left part shows the data for 1 < kL/ (2tt) < 20, and the right part shows the 
data for 20 < kL/(2n) < 64 by changing the x and y axis scaling. In all plots, the solid line is 
the average value for the transfer function. This average was obtained over 161 frames in the 
saturated state and shells of constant |fc|. The upper (lower) dashed line that matches color 
with the solid line correspond to the transfer function plus (minus) one temporal standard 
deviation. From top to bottom, the plots show S (red) and A (black), T bb , T divv , and T bv . 
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Fig. 17. — Kinetic Fourier transfer functions versus a dimensionless wave number (L is 
the length of the smallest dimension in the box) for NZ128. Each plot is displayed in two 
components; the left part shows the data for 1 < kL/ (27r) < 20, and the right part shows the 
data for 20 < kL/ (2tt) < 64 by changing the x and y axis scaling. In all plots, the solid line is 
the average value for the transfer function. This average was obtained over 161 frames in the 
saturated state and shells of constant |fc|. The upper (lower) dashed line that matches color 
with the solid line correspond to the transfer function plus (minus) one temporal standard 



deviation. From top to bottom, the plots show T v 



T vb , T press (red) and T comp (black), and 
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Fig. 18. — Numerical dissipation quantities plotted against a dimensionless wave number (L 
is the length of the smallest dimension in the box). These plots correspond to data from 
NZ128. The upper left plot shows the dissipation rate of kinetic energy (green) and magnetic 
energy (red) in Fourier space. The upper right plot shows the ratio of these two dissipation 
rates. The lower left plot shows the effective numerical viscosity (green) and resistivity (red). 
The lower right plot shows the ratio of the viscosity to resistivity (i.e., the effective Prandtl 
number). In all plots, the solid line is the average value for the quantity of interest. For 
.Dkin and D mag , this average was obtained from averaging over shells of constant |fe| and 
over 161 frames in the saturated state. The averaged viscosity and resistivity values were 
calculated as described in the text. The upper and lower dashed lines correspond to the 
error propagated from one temporal standard deviation. 
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Fig. 19. — Averaged dissipation related quantities as a function of grid resolution. These 
plots correspond to data from the net flux simulations, NZ16, NZ32, NZ64, and NZ128. 
The upper left plot shows the effective viscosity versus x resolution. The dashed line shows 
u c g oc N~ 2 . The upper right plot shows the effective resistivity versus x resolution. Again, 
the dashed line shows i] c g oc N~ 2 . The lower left plot shows the ratio of kinetic to magnetic 
dissipation versus x resolution. The lower right plot shows the effective Prandtl number 
versus x resolution. For each resolution, the data point was obtained from averaging the 
quantity as a function of k over values of k where the error in this quantity is not much larger 
than the mean value. The error bars represent the propagated errors from the temporal 
statistics. 



